engelpath = '../data/DataforMatlab/engel_coefficients.csv';
dataTable = readtable(engelpath);
engel = table2array(dataTable(:,2:52));
engelbet0 = permute(engel(:,+subset),[2,3,1]);
engelbet1 = permute(engel(:,17+subset),[2,3,1]);
engelbet2 = permute(engel(:,34+subset),[2,3,1]); 
b0 = engelbet0 + engelbet1.*log(I_vec) + engelbet2.*log(I_vec).^2;
indT = and(b0>0,b0<1);

eta_m1_Bxi = (engelbet1 + 2*engelbet2.*log(I_vec)).*Bxi_vec./B_vec(subset,:,:);
 
eta_m1_Bxi(indT) = 0; 
sum_eta_m1_Bxi = sum(eta_m1_Bxi);

sigw = 1 - (.075 + Bx_vec.*sum_eta_m1_Bxi)./((1-Bx_vec));  
  
for tt = 1:T
value25 = prctile(I_vec(:,:,tt), 25);
value50 = prctile(I_vec(:,:,tt), 50);
value75 = prctile(I_vec(:,:,tt), 75);

% 
[~, index25] = min(abs(I_vec(:,:,tt) - value25));
[~, index50] = min(abs(I_vec(:,:,tt) - value50));
[~, index75] = min(abs(I_vec(:,:,tt) - value75));

%
sig25(tt) = sigw(1,index25,tt);
sig50(tt) = sigw(1,index50,tt);
sig75(tt) = sigw(1,index75,tt);

sigw(1,1:index25,tt) = sig25(tt);
sigw(1,index75:end,tt) = sig75(tt);
end

sigmaM = sigw;